muldiv 512x512

$$ \gdef\a{\mathtt{a}} \gdef\b{\mathtt{b}} $$

Goal. Compute $$ \mathtt{div512}(\a_0, \a_1, \b_0, \b_1) = \floor{\frac{\a_0 + \a_1 ⋅ 2^{256}}{\b_0 + \b_1 ⋅ 2^{256}}} $$ where all arguments are integers in the range $[0, 2^{256})$.

contract Div512 {
    function div512(uint256 a0, uint256 a1, uint256 b0, uint256 b1)
    internal pure returns (uint256 result)
    {
        if(b1 == 0) {
            return div512x256(a0, a1, b0)
        }
        // ...
    }

Here we can use Knuth's algorithm D, see for example the implementation in OpenZKP.

If $\mathtt{b}_1 = 0$ we can use the previous method, so assume $\mathtt{b}_1 ≥ 1$. In this case the high limbs bound the answer:

$$ \floor{\frac {\mathtt{a}_1} {\mathtt{b}_1 + 1} } ≤ \floor{\frac {\mathtt{a}_0 + \mathtt{a}_1 ⋅ 2^{256}} {\mathtt{b}_0 + \mathtt{b}_1 ⋅ 2^{256}} } ≤ \floor{\frac {\mathtt{a}_1 + 1} {\mathtt{b}_1} } $$

So we can start by subtracting

$$ a' = a - \floor{\frac {\mathtt{a}_1} {\mathtt{b}_1 + 1} } ⋅ b $$

If $a' < b$ we are done.